% This code simulates the model for EIS = 1. Focus on conditional variance
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function SIMDATA = Sim_EIS_1_CV(l_q , l_q_burn , N_sim , N_cond_sim , nu , c_0 , Var_Names)

% INPUT:
% -- l_q: Number of sample periods
% -- l_q_burn: Number of burn-in periods
% -- N_sim: Number of simulations
% -- nu: constant gain parameter
% -- c_0: starting value of log real per capita consumption level
% -- Var_Names: names of variables to be stored
%
% OUTPUT:
% -- SIMDATA: stored simulations as a cell

%% 1. LOAD PARAMETERS
if nu == 0
    error('nu cannot be zero!')
end

I = param_EIS_1(nu);

Constraint = I.tmu_m + 0.5 * (sqrt(1 + I.nu) * (1 + I.nu * (I.lambda - 1)/I.alpha) - I.xi)^2 * I.sigma^2;

if Constraint >= 0
    error('Constraint violated!')
end
%% 2. SIMULATIONS
% additional parameters for simulation
T = l_q;                %Number of sample periods
T_burn = l_q_burn;            %Number of burn-in periods
T_total = T_burn + T;     %Total periods including burn-in periods for simulation
d_0 = c_0 + I.d_c_mean;   %Starting value of log dividend

%% 2.1 Simulate the economy
% Simulate quarterly iid consumption growth, column for each simulation
cg_full = normrnd(I.mu , I.sigma , T_total , N_sim); 

% Simulate quarterly shocks to dividend growth
dg_shock = normrnd(0 , I.sigma_d , T_total , N_sim); 

% Log real consumption level from 1 to T_total
c_full = c_0 + cumsum(cg_full); 

% Log real dividend level from 1 to T_total
d_full = zeros(T_total , N_sim); 

d_full(1 , :) = I.lambda * cg_full(1 , :) + (1 - I.alpha) * d_0 + I.alpha * c_0 + I.alpha * I.mu_dc + dg_shock(1 , :);

for k = 2 : T_total
  d_full(k , :) = I.lambda * cg_full(k , :) + (1 - I.alpha) * d_full(k - 1 , :) + I.alpha * c_full(k - 1 , :) + I.alpha * I.mu_dc + dg_shock(k , :);
end

% Log dividend growth from 1 to T_total
dg_full = [d_full(1 , :) - d_0 ; d_full(2 : end , :) - d_full(1 : end - 1 , :)]; 

% Four-quarter change
cg_yoy_full = (c_full(5 : end , :) - c_full(1 : end - 4 , :))/4;
dg_yoy_full = (d_full(5 : end , :) - d_full(1 : end - 4 , :))/4;

% Construct State Variable \tilde{\mu}
tmu_full = expweightedavg(cg_full , I.nu);
if I.nu == 0
    tmu_full = I.mu * ones(T_total , N_sim);
end

% Calculate Risk-free rate from 1 to T_total
rf_full = -I.tmu_m + tmu_full - 0.5 * I.xi^2 * I.sigma^2; 

Rf_full = exp(rf_full) - 1; 

%% 2.2 Calculate price
P_full = zeros(T_total , N_sim);                 %Price level for the t sample periods

% Max number of terms to calculate explicitly
J = round(log(0.000001)/log(1 - I.alpha)); 

% Calculate the first J terms of A_{c,n}: A_{c,1},...,A_{c,J} as in Equation (D.17)
A_n_inc = ( sqrt(1+I.nu) ... 
          * ( I.nu * (I.lambda-1) * (1 - (1 - I.alpha).^(0 : J - 1))/I.alpha ...
          + (I.lambda - 1) * (1-I.alpha).^(0 : J - 1) + 1 ) ...
          - I.xi ).^2;
      
A_N = cumsum(A_n_inc); % These are the A_{c,n}

% Further incremental terms to A_{c,n} will be approximated by A_n_inc_lim
A_n_inc_lim = ( sqrt(1 + I.nu) * ( I.nu * (I.lambda - 1)/I.alpha + 1) - I.xi )^2;

% Calculate the first J dividend strips' prices at time t. first calculate
% the multiplicative term that is invariant to time on the RHS of Equation (32)
p_cons = ( I.tmu_m * (1 : J) ... % n\tmu_m
       + (1 - (1 - I.alpha).^(1 : J)) * I.mu_dc ...%(1 - (1-\alpha)^n)\mu_{dc}
       + 0.5 * A_N * I.sigma^2 ... %0.5 * A_{c,n} \sigma^2
       + 0.5 * (1 - (1 - I.alpha).^(2 * (1 : J))) * I.sigma_d^2/(1 - (1 - I.alpha)^2) )'; % 0.5 * A_{d,n} \sigma_d^2

% Multiplier that is invariant with respect to simulations or time in
% Equation (D.31) and (D.32)
approx_cons = exp( I.mu_dc + 0.5 * I.sigma_d^2/(1 - (1 - I.alpha)^2) ...
            + (J + 1) * I.tmu_m + 0.5 * A_N(J) * I.sigma^2 + 0.5 * A_n_inc_lim * I.sigma^2 ) ...
            /(1 - exp(I.tmu_m + 0.5 * A_n_inc_lim * I.sigma^2));

% Calculate total price using Equation (D.31)        
for n = 1 : T_total 
    %Calculate p_t^n for maturity up to J
    p_t_n_temp = d_full(n , :) ... %d_t from LHS
  + bsxfun(@times , ( 1 - (1 - I.alpha).^(1 : J) )' , (c_full(n , :) - d_full(n , :) + tmu_full(n,:) * (I.lambda - 1)/I.alpha) ) ...
  + p_cons;

    %Calculate aggregate price
    P_full(n,:) = sum(exp(p_t_n_temp) , 1) ... % First J dividend strip
                + exp(c_full(n , :)) .* exp( (I.lambda - 1) * tmu_full(n,:)/I.alpha ) * approx_cons;    
end

%Calculate trailing 4-quarter sum of dividends
D_ttm_full = movsum(exp(d_full) , 4 , 'Endpoints' , 'discard');

%% 2.3 Calculate other quantities
% P/D ratio
PD_ttm_full = P_full(4 : end , :) ./ D_ttm_full; 
pd_ttm_full = log(PD_ttm_full);

PD_full = P_full./exp(d_full);
pd_full = log(PD_full);

% Objective realized returns
ret_full = log( P_full(2 : end , :) + exp(d_full(2 : end , :)) ) ... % log(P_{t+1} + D_{t+1})
         - log( P_full(1 : end - 1 , :) ); % - log(P_t)
     
ret_MS_full = movsum(ret_full , 4 , 'Endpoints' , 'discard');     
     
Ret_full = exp(ret_full) - 1;

EXRET_full = exp(ret_full) - exp(rf_full(1 : end - 1 , :)); % R_{t+1} - R_{f,t}

exret_full = ret_full - rf_full(1:end-1,:); %r_{t+1} - r_{f,t}

% Experienced dividend growth
tmud_full = expweightedavg(dg_full , I.nu);

% Experience returns
tmur_1_full = expweightedavg(ret_full , I.nu);   %2 to T_total, use experienced log returns(in the paper)
tmur_2_full = expweightedavg(exret_full , I.nu); %2 to T_total, use experienced log excess returns(in case)

%% 3. SIMULATED DATA

%% 3.1 Drop burn-in periods
% All data are re-labled from 1 to T as realizations. 
% Be careful with risk-free rate

% Endowment Process
cg = cg_full(end - T + 1 : end , :);
dg = dg_full(end - T + 1 : end , :);
c  =  c_full(end - T + 1 : end , :);
d  =  d_full(end - T + 1 : end , :);
cg_yoy = cg_yoy_full(end - T + 1 : end , :);
dg_yoy = dg_yoy_full(end - T + 1 : end , :);

% State variable proxies
tmu = tmu_full(end - T + 1 : end , :);
tmud = tmud_full(end - T + 1 : end , :);
tmur_1 = tmur_1_full(end - T + 1 : end , :);
tmur_2 = tmur_2_full(end - T + 1 : end , :);
tmur = tmur_1; % For now we use the experienced log returns specification

% Price Quantity
rf = rf_full(end - T + 1 : end , :);
Rf = Rf_full(end - T + 1 : end , :);

P = P_full(end - T + 1 : end , :);
D_ttm = D_ttm_full(end - T + 1 : end , :);
pd_ttm = pd_ttm_full(end - T + 1 : end , :);
pd = pd_full(end - T + 1 : end , :);

ret = ret_full(end - T + 1 : end , :);
ret_MS = ret_MS_full(end - T + 1 : end , :);
Ret = Ret_full(end - T + 1 : end , :);
EXRET = EXRET_full(end - T + 1 : end , :);
exret = exret_full(end - T + 1 : end , :);

%% 4. CONDITIONAL VARIANCE

condVar = nan(T , N_sim);

for i = 1 : N_sim
    for j = 1 : T
    % Simulate next-period quarterly consumption growth, 1 * N 
    cg_next_temp = normrnd(I.mu , I.sigma , 1 , N_cond_sim); 

    % Simulate next-period quarterly shocks to dividend growth, 1 * N
    dg_shock_next_temp = normrnd(0 , I.sigma_d , 1 , N_cond_sim); 

    % Next-period log real consumption level, 1 * N 
    c_next_temp = c(j , i) + cg_next_temp; 

    % Next-period log real dividend level, 1 * N
    d_next_temp = I.lambda * cg_next_temp + (1 - I.alpha) * d(j , i) + I.alpha * c(j , i) + I.alpha * I.mu_dc + dg_shock_next_temp;
     
    % Next-period state variable, 1 * N 
    tmu_next_temp = (1 - I.nu) * tmu(j , i) + I.nu * cg_next_temp;
    
    % Next period price, 1 * N 
    p_t_n_next_temp = d_next_temp ... %d_t from LHS
  + bsxfun(@times , ( 1 - (1 - I.alpha).^(1 : J) )' , (c_next_temp - d_next_temp + tmu_next_temp * (I.lambda - 1)/I.alpha) ) ...
  + p_cons;

    %Calculate aggregate price, 1 * N 
    P_next_temp = sum(exp(p_t_n_next_temp) , 1) ... % First J dividend strip
                + exp(c_next_temp) .* exp( (I.lambda - 1) * tmu_next_temp/I.alpha ) * approx_cons;  
            
    ret_next_temp = log(P_next_temp + exp(d_next_temp)) ... % log(P_{t+1} + D_{t+1})
                  - log(P(j , i)); % - log(P_t) 
              
    condVar(j , i) = var(ret_next_temp);          
    end
end


%% 5. STORE SIMULATIONS
SIMDATA = cell(length(Var_Names) , 1);
for i = 1 : length(Var_Names)
    SIMDATA{i} = eval(Var_Names{i});
end

end


               